  MODULE TSTORMS_MOD
  implicit none

!=====================================================================
! --- NAMELIST
!=====================================================================

  real :: crit_vort   =   1.6E-4   
  real :: crit_twc    =   0.8   
  real :: crit_thick  =  50.0   
  real :: crit_dist   =   4.0   
  real :: dist_twc    =   5.0
  real :: crit_psl    = 400.0
  real :: dist_psl    =   5.0
  real :: lat_bound_n =  90.0
  real :: lat_bound_s = -90.0

  namelist / nml_tstorms / crit_vort, crit_twc, crit_thick, crit_dist, &
                           dist_twc, crit_psl, dist_psl,               &
                           lat_bound_n, lat_bound_s

!=====================================================================
  contains

!######################################################################

  SUBROUTINE SET_TSTORMS
  implicit none
        READ( *, nml_tstorms )

  end SUBROUTINE SET_TSTORMS

!######################################################################

  SUBROUTINE TSTORMS ( Gwind, Gvort, Gtbar, Gpsl, Gthick,      &
                       Grlon, Grlat, iyear, imon, iday, ihour, &
                       iucy )

!===================================================================
! --- LOCATE TROPICAL STORMS 
!===================================================================

  use TSGPAD_MOD, only : GPAD2,  GPAD1
  use SPLINE_MOD, only : SPLIE2, SPLIE3, FRPRM, SHAPE
  implicit none

!-------------------------------------------------------------------
! --- INPUT ARGUMENTS 
!     Gwind  - wind speed at 850 mb
!     Gvort  - vorticity  at 850 mb
!     Gtbar  - mean temperature for warm core layer
!     Gpsl   - sea level pressure
!     Gthick - thickness of 200 to 1000 mb layer 
!     Grlon  - longitudes
!     Grlat  - latitudes
!     iyear  - year
!     imon   - month
!     iday   - day of month
!     ihour  - hour
!     iucy   - unit for output
!-------------------------------------------------------------------
! --- OUTPUT - file "cyclones" 
!-------------------------------------------------------------------
! --- record # 1 
!     num0   - day
!     imon0  - month
!     iyear  - year
!     number - number of cyclones found
! --- records # 2...number+1 
!     idex, jdex - (i,j) index of cyclone 
!     svort_max  - max vorticity                  
!     swind_max  - max wind              
!      spsl_min  - min sea level pressure                 
!     svort_lon,  svort_lat - longitude & latitude of max vorticity 
!      spsl_lon,   spsl_lat - longitude & latitude of min slp 
!      stwc_lon,   stwc_lat - longitude & latitude of warm core 
!    sthick_lon, sthick_lat - longitude & latitude of max thickness 
!-------------------------------------------------------------------

  real,    intent(in),    dimension(:,:) :: Gwind, Gvort,  Gtbar
  real,    intent(in),    dimension(:,:) :: Gpsl,  Gthick
  real,    intent(in),    dimension(:)   :: Grlon, Grlat
  integer, intent(in)                    :: ihour, iday,  imon, iyear, iucy

!-------------------------------------------------------------------
! --- LOCAL
!-------------------------------------------------------------------

  integer, parameter :: nx   = 16
  integer, parameter :: nx2  = 2*nx
  integer, parameter :: nxp1 = nx + 1

  real,    parameter :: ftol  = 0.01   
  integer, parameter :: nsmax = 10000

  real, dimension(SIZE(Grlon)+nx2) :: rlon
  real, dimension(SIZE(Grlat)+nx2) :: rlat

  real, dimension(SIZE(Gwind,1)+nx2,SIZE(Gwind,2)+nx2) ::  &
        vort, wind, psl,    tbar,    thick,            &
                    psl_dx, tbar_dx, thick_dx,         &
                    psl_dy, tbar_dy, thick_dy 

  real    :: vort_max, wind_max, psl_min, twc_max,   thick_max
  real    :: lon_vort,           lon_psl, lon_twc,   lon_thick
  real    :: lat_vort,           lat_psl, lat_twc,   lat_thick            
  logical ::                              exist_twc, exist_thick

  integer, dimension(nsmax) :: idex, jdex

  real,    dimension(nsmax) ::                                    &
           svort_max, swind_max,  spsl_min, stwc_max, sthick_max, &
           svort_lon,             spsl_lon, stwc_lon, sthick_lon, &
           svort_lat,             spsl_lat, stwc_lat, sthick_lat
  logical, dimension(nsmax) ::                                    &
                                            stwc_is,  sthick_is

  real,    dimension(2)  :: p

  real :: xx, yy, rr, fret

  integer :: ierr_pos, ierr_mag

  integer :: i, im, ip, ix, ixp3, ixp6
  integer :: j, jm, jp, jx, jxp3, jxp6, jxp6h 
  integer :: number, iter
  
!====== Colin stuff

  integer :: ii,jj,psmin_i,psmin_j,psl_crit_dist,starti,startj
  real :: psmn
  
  integer :: kk,ll,tbrmx_i,tbrmx_j, twc_crit_dist,mm, twc_off_dist
  real :: tbar_max, twc_n, twc_s, twc_e, twc_w
  real :: twc_n_temp, twc_s_temp, twc_e_temp, twc_w_temp
  
  integer :: nn,oo,tckmx_i,tckmx_j, thick_crit_dist,pp, thick_off_dist
  real :: tck_max, thick_n, thick_s, thick_e, thick_w
  real :: thick_n_temp, thick_s_temp, thick_e_temp, thick_w_temp
  
  real :: resolution, eps, dist
  integer :: degree_thresh, thresh_grid_dist, center_i, center_j
  integer :: iii, jjj, count
  integer, allocatable :: matrix(:,:), matrix2(:,:), ringmtx(:,:)

!===================================================================

  ix    = SIZE( Gwind, 1 )
  jx    = SIZE( Gwind, 2 )
  ixp3  = ix + nx
  jxp3  = jx + nx
  ixp6  = ix + nx2
  jxp6  = jx + nx2
  jxp6h = jxp6 / 2

  number = 0

!-------------------------------------------------------------------
! --- SETUP
!-------------------------------------------------------------------

  CALL GPAD2( Gwind,  wind  )    !  Wind speed at 850 mb 
  CALL GPAD2( Gvort,  vort  )    !  Vorticity  at 850 mb 
  CALL GPAD2( Gtbar,  tbar  )    !  Mean temp for warm core layer
  CALL GPAD2( Gpsl,   psl   )    !  Sea level pressure.
  CALL GPAD2( Gthick, thick )    !  Thickness of 200 to 1000 mb layer 
  CALL GPAD1( Grlon,  rlon  )    !  Longitudes
  CALL GPAD1( Grlat,  rlat  )    !  Latitudes  
 
! --- change sign of vorticity in southern hemisphere
  do j=1,jxp6
  if( rlat(j) < 0.0 ) then
    vort(:,j) = -1.0 * vort(:,j) 
  end if
  end do

! --- change sign of temperature & thickness
!   tbar(:,:) = -1.0 *  tbar(:,:)
!  thick(:,:) = -1.0 * thick(:,:)

!-------------------------------------------------------------------
! --- INITIALIZE SPLINES
!-------------------------------------------------------------------
 
!  CALL SPLIE2( rlon, rlat, psl,   psl_dy   )
!  CALL SPLIE3( rlon, rlat, psl,   psl_dx   )
 
!  CALL SPLIE2( rlon, rlat, tbar,  tbar_dy  )
!  CALL SPLIE3( rlon, rlat, tbar,  tbar_dx  )
 
!  CALL SPLIE2( rlon, rlat, thick, thick_dy )
!  CALL SPLIE3( rlon, rlat, thick, thick_dx )
  
!-------------------------------------------------------------------
! --- CREATE THRESHOLD MATRICES (COLIN)
!-------------------------------------------------------------------

  resolution = 0.25
  eps = 10e-5

  degree_thresh = 2
  thresh_grid_dist = degree_thresh/resolution

  center_i = thresh_grid_dist+1
  center_j = thresh_grid_dist+1

  allocate ( matrix(2*thresh_grid_dist+1,2*thresh_grid_dist+1) )
  
  do iii = 1,2*thresh_grid_dist+1
    do jjj = 1,2*thresh_grid_dist+1
      dist = ( (iii-center_i)**2 + (jjj-center_j)**2 )
      dist = dist ** (0.5)
      if (dist .LE. (thresh_grid_dist + eps)) then
        matrix(iii,jjj) = 1
      else
        matrix(iii,jjj) = 0
      end if
    end do
  end do
  
!  do iii=1,2*thresh_grid_dist+1
!    print *,( matrix(iii,jjj), jjj=1,2*thresh_grid_dist+1 )
!  enddo
  
!-------------------------------------------------------------------
! --- CREATE DISTANCE MATRICES (COLIN)
!-------------------------------------------------------------------

  degree_thresh = 8
  thresh_grid_dist = degree_thresh/resolution

  center_i = thresh_grid_dist+1
  center_j = thresh_grid_dist+1

  allocate ( matrix2(2*thresh_grid_dist+1,2*thresh_grid_dist+1) )
  
  do iii = 1,2*thresh_grid_dist+1
    do jjj = 1,2*thresh_grid_dist+1
      dist = ( (iii-center_i)**2 + (jjj-center_j)**2 )
      dist = dist ** (0.5)
      if (dist .LE. (thresh_grid_dist + eps)) then
        matrix2(iii,jjj) = 1
      else
        matrix2(iii,jjj) = 0
      end if
    end do
  end do
  

  
!-------------------------------------------------------------------
! --- CREATE RING MATRIX (COLIN)
!-------------------------------------------------------------------

  allocate ( ringmtx(2*thresh_grid_dist+1,2*thresh_grid_dist+1) )
  ringmtx = matrix2
  
  
    do iii = 1,2*thresh_grid_dist+1
        do jjj = 1,2*thresh_grid_dist+1
            count = 1
            if (matrix2(iii,jjj) .EQ. 0) then
                count = count+1
            else
                ringmtx(iii,jjj) = 2
                ringmtx(iii,2*thresh_grid_dist+1-jjj+1) = 2
                exit
            end if
        end do
    end do
    
    do jjj = 1,2*thresh_grid_dist+1
        do iii = 1,2*thresh_grid_dist+1
            count = 1
            if (matrix2(iii,jjj) .EQ. 0) then
                count = count+1
            else if (matrix2(iii,jjj) .EQ. 1) then
                ringmtx(iii,jjj) = 2
                ringmtx(2*thresh_grid_dist+1-iii+1,jjj) = 2
                exit
            else
            end if
        end do
    end do

!  do iii=1,2*thresh_grid_dist+1
!    print *,( ringmtx(iii,jjj), jjj=1,2*thresh_grid_dist+1 )
!  enddo
  
!===================================================================
! -- LOOP OVER GRID & LOOK FOR STORMS
!===================================================================

  do j = nxp1,jxp3
         if( ( rlat(j) > lat_bound_n ) .or. &
             ( rlat(j) < lat_bound_s ) ) CYCLE
  do i = nxp1,ixp3
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz

! Changed from 4 -> 8
    im = i - 16
    ip = i + 16
    jm = j - 16
    jp = j + 16

!-------------------------------------------------------------------
! --- STEP 1: CHECK FOR VORTICITY MAX
!-------------------------------------------------------------------

   vort_max = MAXVAL( vort(im:ip,jm:jp) )
   wind_max = MAXVAL( wind(im:ip,jm:jp) )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  if( ( vort(i,j) /= vort_max  )  .or. &
      ( vort(i,j) <  crit_vort ) ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

  lon_vort = rlon(i)
  lat_vort = rlat(j)

!-------------------------------------------------------------------
! --- STEP 2: LOCATE LOCAL SEA LEVEL PRESSURE MIN 
!-------------------------------------------------------------------

  ierr_pos  = 0

!  p(1) = lon_vort
!  p(2) = lat_vort

!  CALL FRPRM( rlon, rlat, psl,  psl_dy, psl_dx, &
!              p,    ftol, iter, fret,   ierr_pos )
  
! Set large number for initial psmn
  psmn = 1E7
  psmin_i = 1
  psmin_j = 1
 
! Crit dist in gridboxes currently
! 4 for 2 x 2 data
  psl_crit_dist = 8

  starti = i-psl_crit_dist
  startj = j-psl_crit_dist
  
! Loop over gridboxes within critdist to see where psl_min is
  do ii = i-psl_crit_dist,i+psl_crit_dist
  do jj = j-psl_crit_dist,j+psl_crit_dist
! Makes sure it's within circular domain defined by matrix
!  print *,ii-starti+1,' ',jj-startj+1
  if (matrix(ii-starti+1,jj-startj+1) .gt. 0) then
      if (psl(ii,jj).lt.psmn) then
        psmn = psl(ii,jj)
        psmin_i = ii
        psmin_j = jj
      else
      end if
  end if
  end do
  end do
  
  lon_psl = rlon(psmin_i)
  lat_psl = rlat(psmin_j)
  
  if (psmin_i == i-psl_crit_dist) then
  	ierr_pos = 1
  else if (psmin_i == i+psl_crit_dist) then
    ierr_pos = 1
  else if (psmin_j == j-psl_crit_dist) then
    ierr_pos = 1
  else if (psmin_j == j+psl_crit_dist) then
    ierr_pos = 1
  else
    ierr_pos = 0
  end if
  
  psl_min = psmn

!  psl_min = fret
!  lon_psl = p(1)
!  lat_psl = p(2)
!  xx      = lon_psl - lon_vort
!  yy      = lat_psl - lat_vort
!  rr      = xx*xx + yy*yy

!  if( rr >= crit_dist  ) ierr_pos = 1
!  if( psl_min < 500.e2 ) ierr_pos = 1

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

! --- magnitude

  ierr_mag = 0

!  CALL SHAPE( rlon, rlat, psl,  psl_dy, psl_dx,    &
!              p,    ftol, iter, fret,   ierr_mag,  &
!              crit_psl, dist_psl )

! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_mag == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


!-------------------------------------------------------------------
! --- STEP 3: CHECK FOR PRESENCE OF A WARM CORE
!-------------------------------------------------------------------

! --- location

  ierr_pos = 0

!  CALL FRPRM( rlon, rlat, tbar, tbar_dy, tbar_dx, &
!              p,    ftol, iter, fret,    ierr_pos )

!  if( ierr_pos == 0 ) then
!    twc_max = -fret
!    lon_twc = p(1)
!    lat_twc = p(2)
!    xx      = lon_twc - lon_psl
!    yy      = lat_twc - lat_psl
!    rr      = xx*xx + yy*yy
!    if( rr >= crit_dist ) ierr_pos = 1
!  else
!    twc_max = 0.0
!    lon_twc = 0.0
!    lat_twc = 0.0
!  endif
  
  ! 4 for 2 x 2
  twc_crit_dist = 8
  
  tbar_max = 0
  tbrmx_i = 1
  tbrmx_j = 1
  
  starti = i-twc_crit_dist
  startj = j-twc_crit_dist
  
  do kk = i-twc_crit_dist,i+twc_crit_dist
  do ll = j-twc_crit_dist,j+twc_crit_dist
  if (matrix(kk-starti+1,ll-startj+1) .gt. 0) then
      if (tbar(kk,ll).gt.tbar_max) then
        tbar_max = tbar(kk,ll)
        tbrmx_i = kk
        tbrmx_j = ll
      else
      end if
  end if
  end do
  end do

  
  lon_twc = rlon(tbrmx_i)
  lat_twc = rlat(tbrmx_j)
  
  if (tbrmx_i == i-twc_crit_dist) then
  	ierr_pos = 1
  else if (tbrmx_i == i+twc_crit_dist) then
    ierr_pos = 1
  else if (tbrmx_j == j-twc_crit_dist) then
    ierr_pos = 1
  else if (tbrmx_j == j+twc_crit_dist) then
    ierr_pos = 1
  else
    ierr_pos = 0
  end if
  
  twc_max = tbar_max
  

  exist_twc = (ierr_pos == 0)



! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
!   if ( ierr_pos == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

! --- magnitude

  if ( exist_twc ) then

  ierr_mag = 0

!  CALL SHAPE( rlon, rlat, tbar, tbar_dy, tbar_dx,   &
!              p,    ftol, iter, fret,    ierr_mag,  &
!              crit_twc, dist_twc )
! 

! Critical distance in gridboxes which is acceptable for warm core falloff
! 5 for 2 x 2
  twc_off_dist = 20
  
! This is a first guess at warm core falloff by looking at max acceptable radius

  twc_n = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i - twc_off_dist, tbrmx_j)
  twc_s = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i + twc_off_dist, tbrmx_j)
  twc_e = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i, tbrmx_j - twc_off_dist)
  twc_w = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i, tbrmx_j + twc_off_dist)
  
! Loop over radii inside of max acceptable radius to see if there are any other values
! which are deemed acceptable falloff

  do mm = 1,twc_off_dist
  twc_n_temp = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i - mm, tbrmx_j)
  twc_s_temp = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i + mm, tbrmx_j)
  twc_e_temp = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i, tbrmx_j - mm)
  twc_w_temp = tbar(tbrmx_i,tbrmx_j) - tbar(tbrmx_i, tbrmx_j + mm)
  if (twc_n_temp.gt.twc_n) then
  	twc_n = twc_n_temp
  end if
  if (twc_s_temp.gt.twc_s) then
  	twc_s = twc_s_temp
  end if
  if (twc_e_temp.gt.twc_e) then
  	twc_e = twc_e_temp
  end if
  if (twc_w_temp.gt.twc_w) then
  	twc_w = twc_w_temp
  end if
  end do
  
! If > -0.5 (if < 0.5), then reject cyclone
  
  if (twc_n.lt.0.8) then
  	ierr_mag = 1
  else if (twc_s.lt.0.8) then
    ierr_mag = 1
  else if (twc_e.lt.0.8) then
    ierr_mag = 1
  else if (twc_w.lt.0.8) then
    ierr_mag = 1
  else
    ierr_mag = 0
  end if
  
  if (tbar(tbrmx_i,tbrmx_j) .lt. 244) then
    ierr_mag = 1
  end if
  
  print *, lon_psl, lat_psl, '     ', lon_twc, lat_twc
!  print *, tbar(tbrmx_i,tbrmx_j)
!  print *, twc_n,' ',twc_s,' ',twc_e,' ',twc_w
!  print *, exist_thick, ' ',exist_twc, ' ',ierr_mag
!  print *, ' '

  exist_twc = exist_twc .and. ( ierr_mag == 0 )
 
  endif
  
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
   if ( ierr_mag == 1 ) CYCLE
! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!!-------------------------------------------------------------------
!! --- STEP 4: CHECK FOR THICKNESS MAX
!!-------------------------------------------------------------------
!
!! --- location
!
  ierr_pos = 0
!
!!  CALL FRPRM( rlon, rlat, thick, thick_dy, thick_dx,  &
!!              p,    ftol, iter,  fret,     ierr_pos )
!
!!  if( ierr_pos == 0 ) then
!!    thick_max = -fret
!!    lon_thick = p(1)
!!    lat_thick = p(2)
!!    xx        = lon_thick - lon_psl
!!    yy        = lat_thick - lat_psl
!!    rr        = xx*xx + yy*yy
!!    if( rr >= crit_dist ) ierr_pos = 1
!!  else
!    thick_max = 0.0
!    lon_thick = 0.0
!    lat_thick = 0.0
!!  endif
!
!  thick_crit_dist = 8
!  
!  tck_max = 0
!  tckmx_i = 1
!  tckmx_j = 1
!  
!  do nn = i-thick_crit_dist,i+thick_crit_dist
!  do oo = j-thick_crit_dist,j+thick_crit_dist
!  if (thick(nn,oo).lt.tck_max) then
!  	tck_max = thick(nn,oo)
!  	tckmx_i = nn
!  	tckmx_j = oo
!  else
!  end if
!  end do
!  end do
!  
!  lon_thick = rlon(tckmx_i)
!  lat_thick = rlat(tckmx_j)
!  
!  if (tckmx_i == i-thick_crit_dist) then
!  	ierr_pos = 1
!  else if (tckmx_i == i+thick_crit_dist) then
!    ierr_pos = 1
!  else if (tckmx_j == j-thick_crit_dist) then
!    ierr_pos = 1
!  else if (tckmx_j == j+thick_crit_dist) then
!    ierr_pos = 1
!  else
!    ierr_pos = 0
!  end if
!  
!  thick_max = tck_max
!
  exist_thick = (ierr_pos == 0)
!
!! --- magnitude
!
  if ( exist_thick ) then
!
  ierr_mag = 0
!
!!  CALL SHAPE( rlon, rlat, thick, thick_dy, thick_dx,  &
!!              p,    ftol, iter,  fret,     ierr_mag,  &
!!              crit_thick, dist_twc )
!
!  thick_off_dist = 32
!  
!! This is a first guess at warm core falloff by looking at max acceptable radius
!
!  thick_n = thick(tckmx_i,tckmx_j) - thick(tckmx_i - thick_off_dist, tckmx_j)
!  thick_s = thick(tckmx_i,tckmx_j) - thick(tckmx_i + thick_off_dist, tckmx_j)
!  thick_e = thick(tckmx_i,tckmx_j) - thick(tckmx_i, tckmx_j - thick_off_dist)
!  thick_w = thick(tckmx_i,tckmx_j) - thick(tckmx_i, tckmx_j + thick_off_dist)
!  
!! Loop over radii inside of max acceptable radius to see if there are any other values
!! which are deemed acceptable falloff
!
!  do pp = 1,thick_off_dist
!  thick_n_temp = thick(tckmx_i,tckmx_j) - thick(tckmx_i - pp, tckmx_j)
!  thick_s_temp = thick(tckmx_i,tckmx_j) - thick(tckmx_i + pp, tckmx_j)
!  thick_e_temp = thick(tckmx_i,tckmx_j) - thick(tckmx_i, tckmx_j - pp)
!  thick_w_temp = thick(tckmx_i,tckmx_j) - thick(tckmx_i, tckmx_j + pp)
!  if (thick_n_temp.lt.thick_n) then
!  	thick_n = thick_n_temp
!  end if
!  if (thick_s_temp.lt.thick_s) then
!  	thick_s = thick_s_temp
!  end if
!  if (thick_e_temp.lt.thick_e) then
!  	thick_e = thick_e_temp
!  end if
!  if (thick_w_temp.lt.thick_w) then
!  	thick_w = thick_w_temp
!  end if
!  end do
!  
!! If > -0.5 (if < 0.5), then reject cyclone
!  
!  if (thick_n.gt.-50) then
!  	ierr_mag = 1
!  else if (thick_s.gt.-50) then
!    ierr_mag = 1
!  else if (thick_e.gt.-50) then
!    ierr_mag = 1
!  else if (thick_w.gt.-50) then
!    ierr_mag = 1
!  else
!    ierr_mag = 0
!  end if
!
  exist_thick = exist_thick .and. ( ierr_mag == 0 )
! 
  endif
!  
!! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
!   if ( ierr_mag == 1 ) CYCLE
!! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

!-------------------------------------------------------------------
! --- WE HAVE A TROPICAL STORM. SAVE INFO ABOUT STORM
!-------------------------------------------------------------------

    number = number + 1

  if( number > nsmax ) then
    PRINT *, '***************************************'
    PRINT *, '  GOT TOO MANY STORMS - INCREASE nsmax '
    PRINT *, '***************************************'
    STOP
  endif

        idex(number) = i - nx
        jdex(number) = j - nx

   svort_max(number) =  vort_max
   svort_lon(number) =  lon_vort
   svort_lat(number) =  lat_vort

   swind_max(number) =  wind_max

    spsl_min(number) =  psl_min
    spsl_lon(number) =  lon_psl
    spsl_lat(number) =  lat_psl

    stwc_max(number) =  twc_max
    stwc_lon(number) =  lon_twc
    stwc_lat(number) =  lat_twc
    stwc_is (number) =  exist_twc

  sthick_max(number) =  thick_max
  sthick_lon(number) =  lon_thick
  sthick_lat(number) =  lat_thick
  sthick_is (number) =  exist_thick

!  WRITE(*,*) i-nx, j-nx
! zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
  end do
  end do

!===================================================================
! ---  OUTPUT
!===================================================================

  WRITE(iucy,*) iday,  imon, iyear, number, ihour
  WRITE(   *,*) iyear, imon, iday,  ihour,  number

  if( number == 0 ) RETURN

! ---  sthick_is(1:number) = .true.

  do i = 1,number
    WRITE(iucy,*) idex(i),       jdex(i),              & 
              spsl_lon(i),   spsl_lat(i),              &
              swind_max(i), svort_max(i), spsl_min(i), &
                stwc_is(i), sthick_is(i)  
  end do
  
  deallocate(matrix)
  deallocate(matrix2)
  deallocate(ringmtx)
  
!===================================================================
  end SUBROUTINE TSTORMS

!######################################################################
  end MODULE TSTORMS_MOD

